"""
Phase 13: Reconciliation Tests — Do Phase 12 Results Contradict Earlier Papers?
================================================================================
Phase 12 found within EMU:
  - Trade balance dominates (-233***), income balance null
  - Investment dominates S-I (+171***), savings null

Earlier papers found (in general panels):
  - Papers 9A/13: Demographics work through income/non-trade residual, trade null
  - Paper 1: Neither S nor I alone significant (GE co-determination)

RECONCILIATION APPROACH:
  Test 1 (Channel Switching): Run CA decomposition on the FULL panel with
          EMU/CFA interaction terms. If Z₁×EMU is significant on trade but
          Z₁ alone is significant on income, the channels are regime-contingent
          (not contradictory).

  Test 2 (S-I Regime Contingency): Run S vs I on the full panel with union
          interactions. If Z₁×EMU flips the active margin from neither→investment,
          then monetary union changes the adjustment mechanism.

Both tests use the trilemma panel (140 countries) to ensure comparability
with Phase 12 subsamples.
"""

import pandas as pd
import numpy as np
from pathlib import Path
import sys
import json
import warnings
warnings.filterwarnings('ignore')

PROJECT_DIR = Path(__file__).resolve().parent.parent
ROOT_DIR = PROJECT_DIR.parent
sys.path.insert(0, str(ROOT_DIR / "multilateral" / "src"))
from model import PanelGLS

DATA = PROJECT_DIR / "data" / "processed"
OUT_TABLES = PROJECT_DIR / "output" / "tables"
OUT_TABLES.mkdir(parents=True, exist_ok=True)


# ── Union definitions (same as phase 12) ──────────────────────────────

EUROZONE_JOIN = {
    'AUT': 1999, 'BEL': 1999, 'FIN': 1999, 'FRA': 1999, 'DEU': 1999,
    'IRL': 1999, 'ITA': 1999, 'LUX': 1999, 'NLD': 1999, 'PRT': 1999,
    'ESP': 1999, 'GRC': 2001, 'SVN': 2007, 'CYP': 2008, 'MLT': 2008,
    'SVK': 2009, 'EST': 2011, 'LVA': 2014, 'LTU': 2015,
}

CFA_JOIN = {
    'BEN': 1960, 'BFA': 1960, 'CIV': 1960, 'MLI': 1984,
    'NER': 1960, 'SEN': 1960, 'TGO': 1960, 'GNB': 1997,
    'CMR': 1960, 'CAF': 1960, 'TCD': 1960, 'COG': 1960,
    'GNQ': 1985, 'GAB': 1960,
}


def stars(p):
    if p < 0.01: return '***'
    if p < 0.05: return '**'
    if p < 0.1: return '*'
    return ''


def fmt(val, se, p):
    return f"{val:.3f}{stars(p)}", f"({se:.3f})"


def run_gls(df, y_var, x_vars, label):
    """Run PanelGLS and return results dict."""
    cols = [y_var] + x_vars + ['iso3', 'year']
    sub = df[cols].dropna()
    if len(sub) < 30:
        print(f"  {label}: insufficient obs ({len(sub)}), skipping")
        return None
    gls = PanelGLS()
    y = sub[y_var].values
    X = sub[x_vars].values
    try:
        gls.fit(y, X, sub['iso3'].values, sub['year'].values)
    except Exception as e:
        print(f"  {label}: GLS failed ({e}), skipping")
        return None
    result = {
        'model': label, 'dep_var': y_var,
        'n_obs': gls.n_obs, 'n_countries': gls.n_countries,
        'r_squared': gls.r_squared,
    }
    for i, var in enumerate(x_vars):
        result[f'{var}_coef'] = gls.beta[i]
        result[f'{var}_se'] = gls.se[i]
        result[f'{var}_p'] = gls.pvalues[i]
    return result


def write_table(rows, headers, filepath, title, notes=None):
    """Write markdown table."""
    with open(filepath, 'w') as f:
        f.write(f"# {title}\n\n")
        f.write('| ' + ' | '.join(headers) + ' |\n')
        f.write('|' + '|'.join(['---'] * len(headers)) + '|\n')
        for row in rows:
            f.write('| ' + ' | '.join(str(c) for c in row) + ' |\n')
        if notes:
            f.write(f"\n{notes}\n")
    print(f"  Wrote: {filepath.name}")


# ── Data loading ──────────────────────────────────────────────────────

def download_ca_components(tri):
    """Download trade balance, income balance, secondary income from WDI."""
    import urllib.request

    indicators = {
        'BN.GSR.GNFS.CD': 'trade_balance_usd',
        'BN.GSR.FCTY.CD': 'primary_income_usd',
        'BN.TRF.CURR.CD': 'secondary_income_usd',
    }

    all_dfs = []
    for code, name in indicators.items():
        json_url = f"https://api.worldbank.org/v2/country/all/indicator/{code}?format=json&per_page=20000&date=1995:2024"
        try:
            req = urllib.request.Request(json_url)
            req.add_header('User-Agent', 'Mozilla/5.0')
            with urllib.request.urlopen(req, timeout=30) as resp:
                data = json.loads(resp.read().decode())
                if len(data) < 2 or data[1] is None:
                    print(f"  WDI {code}: no data returned")
                    continue
                records = []
                for item in data[1]:
                    if item['value'] is not None:
                        records.append({
                            'iso3': item['countryiso3code'],
                            'year': int(item['date']),
                            name: float(item['value'])
                        })
                df = pd.DataFrame(records)
                all_dfs.append(df)
                print(f"  WDI {code} ({name}): {len(df)} obs")
        except Exception as e:
            print(f"  WDI {code} failed: {e}")

    if not all_dfs:
        return None

    result = all_dfs[0]
    for df in all_dfs[1:]:
        result = result.merge(df, on=['iso3', 'year'], how='outer')

    gdp = tri[['iso3', 'year', 'ngdp_usd']].dropna()
    result = result.merge(gdp, on=['iso3', 'year'], how='inner')

    for usd_col in ['trade_balance_usd', 'primary_income_usd', 'secondary_income_usd']:
        gdp_col = usd_col.replace('_usd', '_gdp')
        if usd_col in result.columns:
            result[gdp_col] = (result[usd_col] / (result['ngdp_usd'] * 1e9)) * 100

    result = result.drop(columns=['ngdp_usd'] + [c for c in result.columns if c.endswith('_usd')], errors='ignore')
    return result


def load_data():
    """Load trilemma panel with CA components and union flags."""
    print("Loading trilemma panel...")
    tri = pd.read_csv(DATA / "trilemma_panel.csv")

    # Download CA components
    print("Downloading CA components from WDI...")
    ca_comp = download_ca_components(tri)
    if ca_comp is not None:
        tri = tri.merge(ca_comp, on=['iso3', 'year'], how='left')
        print(f"  CA components merged: {[c for c in ca_comp.columns if c != 'iso3' and c != 'year']}")

    # Build S-I gap
    tri['si_gap'] = tri['gross_savings_gdp'] - tri['gross_investment_gdp']
    tri['non_si_residual'] = tri['ca_gdp'] - tri['si_gap']

    # Union membership flags (time-varying)
    tri['in_emu'] = 0
    for iso, yr in EUROZONE_JOIN.items():
        tri.loc[(tri['iso3'] == iso) & (tri['year'] >= yr), 'in_emu'] = 1

    tri['in_cfa'] = 0
    for iso, yr in CFA_JOIN.items():
        tri.loc[(tri['iso3'] == iso) & (tri['year'] >= yr), 'in_cfa'] = 1

    tri['in_union'] = ((tri['in_emu'] == 1) | (tri['in_cfa'] == 1)).astype(int)

    # Interaction terms
    for z in ['Z_1', 'Z_2', 'Z_3']:
        tri[f'{z}_x_emu'] = tri[z] * tri['in_emu']
        tri[f'{z}_x_cfa'] = tri[z] * tri['in_cfa']
        tri[f'{z}_x_union'] = tri[z] * tri['in_union']

    # Restrict to post-1995 with macro data
    tri = tri[tri['year'] >= 1995].copy()

    n_emu = tri[tri['in_emu'] == 1]['iso3'].nunique()
    n_cfa = tri[tri['in_cfa'] == 1]['iso3'].nunique()
    n_total = tri['iso3'].nunique()
    print(f"  Panel: {len(tri)} obs, {n_total} countries ({n_emu} EMU, {n_cfa} CFA)")

    return tri


# ══════════════════════════════════════════════════════════════════════
# TEST 1: Channel Switching — Trade vs Income with Union Interactions
# ══════════════════════════════════════════════════════════════════════

def test1_channel_switching(tri):
    """
    Run Z → each CA component on FULL panel with Z×EMU and Z×CFA interactions.

    If Z₁ alone is significant on income (replicating Papers 9A/13) and
    Z₁×EMU is significant on trade (replicating Phase 12), the channels
    are regime-contingent: monetary union switches the channel from income
    to trade.
    """
    print("\n" + "=" * 70)
    print("TEST 1: Channel Switching — CA Decomposition with Union Interactions")
    print("=" * 70)

    controls = ['fiscal_bal_gdp', 'nfa_gdp_lag', 'rgdp_growth', 'log_rel_opw', 'kaopen']

    # Dependent variables for decomposition
    dvs = {
        'ca_gdp': 'CA/GDP',
        'trade_balance_gdp': 'Trade Balance/GDP',
        'primary_income_gdp': 'Primary Income/GDP',
        'secondary_income_gdp': 'Secondary Income/GDP',
        'si_gap': 'S-I Gap',
        'non_si_residual': 'Non-S-I Residual',
    }

    results = []

    for dv, dv_label in dvs.items():
        if dv not in tri.columns:
            print(f"  {dv}: not available, skipping")
            continue

        # Model A: Baseline (Z only, no interactions) — replicates earlier papers
        x_base = ['Z_1', 'Z_2', 'Z_3'] + controls
        r_base = run_gls(tri, dv, x_base, f"base_{dv}")
        if r_base:
            r_base['test'] = 'baseline'
            r_base['dv_label'] = dv_label
            results.append(r_base)
            print(f"  Base {dv}: Z₁={r_base['Z_1_coef']:.1f} (p={r_base['Z_1_p']:.4f})")

        # Model B: With EMU interaction
        x_emu = ['Z_1', 'Z_2', 'Z_3', 'Z_1_x_emu', 'Z_2_x_emu', 'Z_3_x_emu', 'in_emu'] + controls
        r_emu = run_gls(tri, dv, x_emu, f"emu_int_{dv}")
        if r_emu:
            r_emu['test'] = 'emu_interaction'
            r_emu['dv_label'] = dv_label
            results.append(r_emu)
            z1_base = r_emu['Z_1_coef']
            z1_int = r_emu['Z_1_x_emu_coef']
            z1_emu_total = z1_base + z1_int
            print(f"  EMU int {dv}: Z₁={z1_base:.1f} (p={r_emu['Z_1_p']:.4f}), "
                  f"Z₁×EMU={z1_int:.1f} (p={r_emu['Z_1_x_emu_p']:.4f}), "
                  f"total EMU={z1_emu_total:.1f}")

        # Model C: With CFA interaction
        x_cfa = ['Z_1', 'Z_2', 'Z_3', 'Z_1_x_cfa', 'Z_2_x_cfa', 'Z_3_x_cfa', 'in_cfa'] + controls
        r_cfa = run_gls(tri, dv, x_cfa, f"cfa_int_{dv}")
        if r_cfa:
            r_cfa['test'] = 'cfa_interaction'
            r_cfa['dv_label'] = dv_label
            results.append(r_cfa)
            z1_base = r_cfa['Z_1_coef']
            z1_int = r_cfa['Z_1_x_cfa_coef']
            print(f"  CFA int {dv}: Z₁={z1_base:.1f} (p={r_cfa['Z_1_p']:.4f}), "
                  f"Z₁×CFA={z1_int:.1f} (p={r_cfa['Z_1_x_cfa_p']:.4f})")

        # Model D: Joint (both EMU and CFA interactions)
        x_joint = (['Z_1', 'Z_2', 'Z_3',
                     'Z_1_x_emu', 'Z_2_x_emu', 'Z_3_x_emu', 'in_emu',
                     'Z_1_x_cfa', 'Z_2_x_cfa', 'Z_3_x_cfa', 'in_cfa']
                    + controls)
        r_joint = run_gls(tri, dv, x_joint, f"joint_{dv}")
        if r_joint:
            r_joint['test'] = 'joint_interaction'
            r_joint['dv_label'] = dv_label
            results.append(r_joint)

    # ── Write Table ──────────────────────────────────────────────────
    if results:
        _write_test1_table(results)
        _write_test1_interpretation(results)

    return results


def _write_test1_table(results):
    """Write the main channel-switching results table."""
    headers = ['DV', 'Spec', 'Z₁ (base)', '', 'Z₁×EMU', '', 'Z₁×CFA', '', 'N', 'R²']
    rows = []

    # Group by dep var
    by_dv = {}
    for r in results:
        key = r['dep_var']
        if key not in by_dv:
            by_dv[key] = {}
        by_dv[key][r['test']] = r

    for dv in ['ca_gdp', 'trade_balance_gdp', 'primary_income_gdp',
               'secondary_income_gdp', 'si_gap', 'non_si_residual']:
        if dv not in by_dv:
            continue

        specs = by_dv[dv]
        dv_label = specs[list(specs.keys())[0]].get('dv_label', dv)

        # Baseline
        if 'baseline' in specs:
            r = specs['baseline']
            c, s = fmt(r['Z_1_coef'], r['Z_1_se'], r['Z_1_p'])
            rows.append([dv_label, 'Base', c, s, '', '', '', '', r['n_obs'], f"{r['r_squared']:.3f}"])

        # EMU interaction
        if 'emu_interaction' in specs:
            r = specs['emu_interaction']
            c1, s1 = fmt(r['Z_1_coef'], r['Z_1_se'], r['Z_1_p'])
            c2, s2 = fmt(r['Z_1_x_emu_coef'], r['Z_1_x_emu_se'], r['Z_1_x_emu_p'])
            rows.append([dv_label, '+EMU', c1, s1, c2, s2, '', '', r['n_obs'], f"{r['r_squared']:.3f}"])

        # CFA interaction
        if 'cfa_interaction' in specs:
            r = specs['cfa_interaction']
            c1, s1 = fmt(r['Z_1_coef'], r['Z_1_se'], r['Z_1_p'])
            c3, s3 = fmt(r['Z_1_x_cfa_coef'], r['Z_1_x_cfa_se'], r['Z_1_x_cfa_p'])
            rows.append([dv_label, '+CFA', c1, s1, '', '', c3, s3, r['n_obs'], f"{r['r_squared']:.3f}"])

        # Joint
        if 'joint_interaction' in specs:
            r = specs['joint_interaction']
            c1, s1 = fmt(r['Z_1_coef'], r['Z_1_se'], r['Z_1_p'])
            c2, s2 = fmt(r['Z_1_x_emu_coef'], r['Z_1_x_emu_se'], r['Z_1_x_emu_p'])
            c3, s3 = fmt(r['Z_1_x_cfa_coef'], r['Z_1_x_cfa_se'], r['Z_1_x_cfa_p'])
            rows.append([dv_label, 'Joint', c1, s1, c2, s2, c3, s3, r['n_obs'], f"{r['r_squared']:.3f}"])

        # Blank separator
        rows.append([''] * 10)

    write_table(rows, headers, OUT_TABLES / "phase13_test1_channel_switching.md",
                "Test 1: Channel Switching — Do Monetary Unions Change the Demographic CA Channel?",
                ("*Notes: PanelGLS with country and year FE. Base = Z₁ alone (replicates earlier papers). "
                 "+EMU/+CFA add Z×union interactions. Joint includes both. "
                 "If Z₁ is significant on income in the base but Z₁×EMU is significant on trade, "
                 "the channels are regime-contingent.*\n"
                 "*Controls: fiscal balance, lagged NFA, growth, relative productivity, KAOPEN.*\n"
                 "*\\*p<0.1, \\*\\*p<0.05, \\*\\*\\*p<0.01*"))


def _write_test1_interpretation(results):
    """Write interpretation of test 1 results."""
    by_dv_test = {}
    for r in results:
        by_dv_test[(r['dep_var'], r['test'])] = r

    lines = ["# Test 1 Interpretation: Channel Switching\n"]

    # Check the key hypotheses
    # H1: Z₁ significant on income (base) — replicates Papers 9A/13
    income_base = by_dv_test.get(('primary_income_gdp', 'baseline'))
    nsi_base = by_dv_test.get(('non_si_residual', 'baseline'))
    trade_base = by_dv_test.get(('trade_balance_gdp', 'baseline'))

    lines.append("## Replication of Earlier Papers (Base Spec, Full Panel)\n")
    if trade_base:
        lines.append(f"- Trade balance: Z₁ = {trade_base['Z_1_coef']:.1f} (p = {trade_base['Z_1_p']:.4f})")
    if income_base:
        lines.append(f"- Primary income: Z₁ = {income_base['Z_1_coef']:.1f} (p = {income_base['Z_1_p']:.4f})")
    if nsi_base:
        lines.append(f"- Non-S-I residual: Z₁ = {nsi_base['Z_1_coef']:.1f} (p = {nsi_base['Z_1_p']:.4f})")

    # H2: Z₁×EMU significant on trade — replicates Phase 12
    trade_emu = by_dv_test.get(('trade_balance_gdp', 'emu_interaction'))
    income_emu = by_dv_test.get(('primary_income_gdp', 'emu_interaction'))

    lines.append("\n## EMU Interaction (Channel Switching Test)\n")
    if trade_emu:
        z1_int = trade_emu.get('Z_1_x_emu_coef', np.nan)
        z1_int_p = trade_emu.get('Z_1_x_emu_p', np.nan)
        z1_base = trade_emu['Z_1_coef']
        lines.append(f"- Trade: Z₁ (non-EMU) = {z1_base:.1f}, Z₁×EMU = {z1_int:.1f} (p = {z1_int_p:.4f})")
        lines.append(f"  → EMU total = {z1_base + z1_int:.1f}")
    if income_emu:
        z1_int = income_emu.get('Z_1_x_emu_coef', np.nan)
        z1_int_p = income_emu.get('Z_1_x_emu_p', np.nan)
        z1_base = income_emu['Z_1_coef']
        lines.append(f"- Income: Z₁ (non-EMU) = {z1_base:.1f}, Z₁×EMU = {z1_int:.1f} (p = {z1_int_p:.4f})")

    # Verdict
    lines.append("\n## Verdict\n")
    channel_switch = False
    if trade_emu and income_base:
        trade_int_sig = trade_emu.get('Z_1_x_emu_p', 1) < 0.1
        income_base_sig = income_base.get('Z_1_p', 1) < 0.1 or (nsi_base and nsi_base.get('Z_1_p', 1) < 0.1)
        if trade_int_sig and income_base_sig:
            channel_switch = True
            lines.append("**CHANNEL SWITCHING CONFIRMED**: Demographics work through income/non-trade in the ")
            lines.append("general panel but Z₁×EMU activates the trade channel within monetary union. ")
            lines.append("The Phase 12 and Papers 9A/13 findings are not contradictory — they are ")
            lines.append("regime-contingent. Monetary union doesn't just amplify the demographic CA ")
            lines.append("effect; it changes the channel through which it operates.")
        elif trade_int_sig:
            lines.append("**PARTIAL**: Z₁×EMU activates the trade channel within EMU, but the income ")
            lines.append("channel is weaker than expected in the base specification. Possible composition ")
            lines.append("effect from using the trilemma panel vs. the extensions/multilateral panel.")
        else:
            lines.append("**NOT CONFIRMED**: Z₁×EMU is not significant on trade balance. The Phase 12 ")
            lines.append("within-EMU trade dominance may be a sample-specific artifact of the within-union ")
            lines.append("demeaning approach rather than a genuine channel switch.")

    filepath = OUT_TABLES / "phase13_test1_interpretation.md"
    with open(filepath, 'w') as f:
        f.write('\n'.join(lines) + '\n')
    print(f"  Wrote: {filepath.name}")

    return channel_switch


# ══════════════════════════════════════════════════════════════════════
# TEST 2: S-I Regime Contingency — Does Union Switch Active Margin?
# ══════════════════════════════════════════════════════════════════════

def test2_si_regime_contingency(tri):
    """
    Run Z → S, Z → I on the full panel with Z×EMU and Z×CFA interactions.

    Paper 1 found neither S nor I alone significant (co-determination).
    Phase 12 found investment dominates within EMU (+171***), savings null.

    If Z₁×EMU is significant on investment but Z₁ alone is null on both
    S and I, then monetary union forces the demographic adjustment onto
    the investment margin specifically.
    """
    print("\n" + "=" * 70)
    print("TEST 2: S-I Regime Contingency — Union Interaction on S vs I")
    print("=" * 70)

    controls = ['fiscal_bal_gdp', 'nfa_gdp_lag', 'rgdp_growth', 'log_rel_opw', 'kaopen']

    dvs = {
        'gross_savings_gdp': 'Gross Savings/GDP',
        'gross_investment_gdp': 'Gross Investment/GDP',
        'si_gap': 'S-I Gap',
    }

    results = []

    for dv, dv_label in dvs.items():
        if dv not in tri.columns:
            print(f"  {dv}: not available, skipping")
            continue

        # Model A: Baseline (no interaction) — replicates Paper 1
        x_base = ['Z_1', 'Z_2', 'Z_3'] + controls
        r_base = run_gls(tri, dv, x_base, f"base_{dv}")
        if r_base:
            r_base['test'] = 'baseline'
            r_base['dv_label'] = dv_label
            results.append(r_base)
            print(f"  Base {dv}: Z₁={r_base['Z_1_coef']:.1f} (p={r_base['Z_1_p']:.4f})")

        # Model B: EMU interaction
        x_emu = ['Z_1', 'Z_2', 'Z_3', 'Z_1_x_emu', 'Z_2_x_emu', 'Z_3_x_emu', 'in_emu'] + controls
        r_emu = run_gls(tri, dv, x_emu, f"emu_int_{dv}")
        if r_emu:
            r_emu['test'] = 'emu_interaction'
            r_emu['dv_label'] = dv_label
            results.append(r_emu)
            z1_base = r_emu['Z_1_coef']
            z1_int = r_emu['Z_1_x_emu_coef']
            print(f"  EMU int {dv}: Z₁={z1_base:.1f} (p={r_emu['Z_1_p']:.4f}), "
                  f"Z₁×EMU={z1_int:.1f} (p={r_emu['Z_1_x_emu_p']:.4f})")

        # Model C: CFA interaction
        x_cfa = ['Z_1', 'Z_2', 'Z_3', 'Z_1_x_cfa', 'Z_2_x_cfa', 'Z_3_x_cfa', 'in_cfa'] + controls
        r_cfa = run_gls(tri, dv, x_cfa, f"cfa_int_{dv}")
        if r_cfa:
            r_cfa['test'] = 'cfa_interaction'
            r_cfa['dv_label'] = dv_label
            results.append(r_cfa)
            z1_base = r_cfa['Z_1_coef']
            z1_int = r_cfa['Z_1_x_cfa_coef']
            print(f"  CFA int {dv}: Z₁={z1_base:.1f} (p={r_cfa['Z_1_p']:.4f}), "
                  f"Z₁×CFA={z1_int:.1f} (p={r_cfa['Z_1_x_cfa_p']:.4f})")

        # Model D: Joint
        x_joint = (['Z_1', 'Z_2', 'Z_3',
                     'Z_1_x_emu', 'Z_2_x_emu', 'Z_3_x_emu', 'in_emu',
                     'Z_1_x_cfa', 'Z_2_x_cfa', 'Z_3_x_cfa', 'in_cfa']
                    + controls)
        r_joint = run_gls(tri, dv, x_joint, f"joint_{dv}")
        if r_joint:
            r_joint['test'] = 'joint_interaction'
            r_joint['dv_label'] = dv_label
            results.append(r_joint)

    # ── Write Table ──────────────────────────────────────────────────
    if results:
        _write_test2_table(results)
        _write_test2_interpretation(results)

    return results


def _write_test2_table(results):
    """Write the S-I regime contingency table."""
    headers = ['DV', 'Spec', 'Z₁ (base)', '', 'Z₁×EMU', '', 'Z₁×CFA', '', 'N', 'R²']
    rows = []

    by_dv = {}
    for r in results:
        key = r['dep_var']
        if key not in by_dv:
            by_dv[key] = {}
        by_dv[key][r['test']] = r

    for dv in ['gross_savings_gdp', 'gross_investment_gdp', 'si_gap']:
        if dv not in by_dv:
            continue

        specs = by_dv[dv]
        dv_label = specs[list(specs.keys())[0]].get('dv_label', dv)

        if 'baseline' in specs:
            r = specs['baseline']
            c, s = fmt(r['Z_1_coef'], r['Z_1_se'], r['Z_1_p'])
            rows.append([dv_label, 'Base', c, s, '', '', '', '', r['n_obs'], f"{r['r_squared']:.3f}"])

        if 'emu_interaction' in specs:
            r = specs['emu_interaction']
            c1, s1 = fmt(r['Z_1_coef'], r['Z_1_se'], r['Z_1_p'])
            c2, s2 = fmt(r['Z_1_x_emu_coef'], r['Z_1_x_emu_se'], r['Z_1_x_emu_p'])
            rows.append([dv_label, '+EMU', c1, s1, c2, s2, '', '', r['n_obs'], f"{r['r_squared']:.3f}"])

        if 'cfa_interaction' in specs:
            r = specs['cfa_interaction']
            c1, s1 = fmt(r['Z_1_coef'], r['Z_1_se'], r['Z_1_p'])
            c3, s3 = fmt(r['Z_1_x_cfa_coef'], r['Z_1_x_cfa_se'], r['Z_1_x_cfa_p'])
            rows.append([dv_label, '+CFA', c1, s1, '', '', c3, s3, r['n_obs'], f"{r['r_squared']:.3f}"])

        if 'joint_interaction' in specs:
            r = specs['joint_interaction']
            c1, s1 = fmt(r['Z_1_coef'], r['Z_1_se'], r['Z_1_p'])
            c2, s2 = fmt(r['Z_1_x_emu_coef'], r['Z_1_x_emu_se'], r['Z_1_x_emu_p'])
            c3, s3 = fmt(r['Z_1_x_cfa_coef'], r['Z_1_x_cfa_se'], r['Z_1_x_cfa_p'])
            rows.append([dv_label, 'Joint', c1, s1, c2, s2, c3, s3, r['n_obs'], f"{r['r_squared']:.3f}"])

        rows.append([''] * 10)

    write_table(rows, headers, OUT_TABLES / "phase13_test2_si_contingency.md",
                "Test 2: S-I Regime Contingency — Does Monetary Union Switch the Active Margin?",
                ("*Notes: PanelGLS with country and year FE. Paper 1 found neither S nor I alone significant "
                 "(co-determination). Phase 12 found investment dominates within EMU. "
                 "If Z₁×EMU activates the investment margin specifically, monetary union forces demographic "
                 "adjustment onto investment.*\n"
                 "*Controls: fiscal balance, lagged NFA, growth, relative productivity, KAOPEN.*\n"
                 "*\\*p<0.1, \\*\\*p<0.05, \\*\\*\\*p<0.01*"))


def _write_test2_interpretation(results):
    """Write interpretation of test 2 results."""
    by_dv_test = {}
    for r in results:
        by_dv_test[(r['dep_var'], r['test'])] = r

    lines = ["# Test 2 Interpretation: S-I Regime Contingency\n"]

    # Check baseline: neither S nor I significant?
    s_base = by_dv_test.get(('gross_savings_gdp', 'baseline'))
    i_base = by_dv_test.get(('gross_investment_gdp', 'baseline'))
    si_base = by_dv_test.get(('si_gap', 'baseline'))

    lines.append("## Replication of Paper 1 (Base Spec, Full Panel)\n")
    if s_base:
        lines.append(f"- Savings: Z₁ = {s_base['Z_1_coef']:.1f} (p = {s_base['Z_1_p']:.4f})")
    if i_base:
        lines.append(f"- Investment: Z₁ = {i_base['Z_1_coef']:.1f} (p = {i_base['Z_1_p']:.4f})")
    if si_base:
        lines.append(f"- S-I gap: Z₁ = {si_base['Z_1_coef']:.1f} (p = {si_base['Z_1_p']:.4f})")

    # Check EMU interaction on investment
    i_emu = by_dv_test.get(('gross_investment_gdp', 'emu_interaction'))
    s_emu = by_dv_test.get(('gross_savings_gdp', 'emu_interaction'))

    lines.append("\n## EMU Interaction (Margin Switching Test)\n")
    if i_emu:
        z1_int = i_emu.get('Z_1_x_emu_coef', np.nan)
        z1_int_p = i_emu.get('Z_1_x_emu_p', np.nan)
        z1_base = i_emu['Z_1_coef']
        lines.append(f"- Investment: Z₁ (non-EMU) = {z1_base:.1f}, Z₁×EMU = {z1_int:.1f} (p = {z1_int_p:.4f})")
        lines.append(f"  → EMU total = {z1_base + z1_int:.1f}")
    if s_emu:
        z1_int = s_emu.get('Z_1_x_emu_coef', np.nan)
        z1_int_p = s_emu.get('Z_1_x_emu_p', np.nan)
        z1_base = s_emu['Z_1_coef']
        lines.append(f"- Savings: Z₁ (non-EMU) = {z1_base:.1f}, Z₁×EMU = {z1_int:.1f} (p = {z1_int_p:.4f})")

    # Verdict
    lines.append("\n## Verdict\n")
    if i_emu and s_emu:
        i_int_sig = i_emu.get('Z_1_x_emu_p', 1) < 0.1
        s_int_sig = s_emu.get('Z_1_x_emu_p', 1) < 0.1

        if i_int_sig and not s_int_sig:
            lines.append("**INVESTMENT MARGIN CONFIRMED**: Z₁×EMU is significant on investment but null on ")
            lines.append("savings. Monetary union forces demographic adjustment onto the investment margin. ")
            lines.append("The Paper 1 'neither alone significant' finding reflects pooled co-determination; ")
            lines.append("within EMU, investment absorbs the full adjustment because exchange rate and ")
            lines.append("interest rate channels are eliminated.")
        elif i_int_sig and s_int_sig:
            lines.append("**BOTH MARGINS ACTIVATE**: Z₁×EMU is significant on both S and I within EMU. ")
            lines.append("Monetary union amplifies the demographic effect on both margins, not just investment.")
        elif not i_int_sig and not s_int_sig:
            lines.append("**NOT CONFIRMED**: Z₁×EMU is null on both S and I. The Phase 12 within-EMU ")
            lines.append("investment dominance may be an artifact of the within-union demeaning approach.")
        else:
            lines.append("**SAVINGS MARGIN**: Unexpected — Z₁×EMU activates savings, not investment.")

    filepath = OUT_TABLES / "phase13_test2_interpretation.md"
    with open(filepath, 'w') as f:
        f.write('\n'.join(lines) + '\n')
    print(f"  Wrote: {filepath.name}")


# ══════════════════════════════════════════════════════════════════════
# SUMMARY
# ══════════════════════════════════════════════════════════════════════

def write_summary(test1_results, test2_results):
    """Write combined summary of both reconciliation tests."""
    print("\n" + "=" * 70)
    print("PHASE 13 SUMMARY: Reconciliation Verdicts")
    print("=" * 70)

    lines = ["# Phase 13: Reconciliation Summary\n"]
    lines.append("## Question: Do Phase 12 findings contradict earlier papers?\n")

    # Collect key results
    t1_by = {(r['dep_var'], r['test']): r for r in test1_results}
    t2_by = {(r['dep_var'], r['test']): r for r in test2_results}

    lines.append("### Test 1: Trade vs Income Channel\n")
    lines.append("| Finding | Source | Z₁ | p |")
    lines.append("|---|---|---|---|")

    # Earlier: income/non-trade significant
    for key_dv, label in [('non_si_residual', 'Non-S-I (base)'),
                           ('primary_income_gdp', 'Income (base)'),
                           ('trade_balance_gdp', 'Trade (base)')]:
        r = t1_by.get((key_dv, 'baseline'))
        if r:
            lines.append(f"| {label} | Full panel | {r['Z_1_coef']:.1f} | {r['Z_1_p']:.4f} |")

    # Phase 12: trade dominates within EMU
    r_trade_emu = t1_by.get(('trade_balance_gdp', 'emu_interaction'))
    if r_trade_emu:
        z_int = r_trade_emu.get('Z_1_x_emu_coef', np.nan)
        p_int = r_trade_emu.get('Z_1_x_emu_p', np.nan)
        lines.append(f"| Trade Z₁×EMU | Interaction | {z_int:.1f} | {p_int:.4f} |")

    r_inc_emu = t1_by.get(('primary_income_gdp', 'emu_interaction'))
    if r_inc_emu:
        z_int = r_inc_emu.get('Z_1_x_emu_coef', np.nan)
        p_int = r_inc_emu.get('Z_1_x_emu_p', np.nan)
        lines.append(f"| Income Z₁×EMU | Interaction | {z_int:.1f} | {p_int:.4f} |")

    lines.append("\n### Test 2: S vs I Margin\n")
    lines.append("| Finding | Source | Z₁ | p |")
    lines.append("|---|---|---|---|")

    for key_dv, label in [('gross_savings_gdp', 'Savings (base)'),
                           ('gross_investment_gdp', 'Investment (base)')]:
        r = t2_by.get((key_dv, 'baseline'))
        if r:
            lines.append(f"| {label} | Full panel | {r['Z_1_coef']:.1f} | {r['Z_1_p']:.4f} |")

    r_i_emu = t2_by.get(('gross_investment_gdp', 'emu_interaction'))
    if r_i_emu:
        z_int = r_i_emu.get('Z_1_x_emu_coef', np.nan)
        p_int = r_i_emu.get('Z_1_x_emu_p', np.nan)
        lines.append(f"| Investment Z₁×EMU | Interaction | {z_int:.1f} | {p_int:.4f} |")

    r_s_emu = t2_by.get(('gross_savings_gdp', 'emu_interaction'))
    if r_s_emu:
        z_int = r_s_emu.get('Z_1_x_emu_coef', np.nan)
        p_int = r_s_emu.get('Z_1_x_emu_p', np.nan)
        lines.append(f"| Savings Z₁×EMU | Interaction | {z_int:.1f} | {p_int:.4f} |")

    filepath = OUT_TABLES / "phase13_summary.md"
    with open(filepath, 'w') as f:
        f.write('\n'.join(lines) + '\n')
    print(f"  Wrote: {filepath.name}")


# ══════════════════════════════════════════════════════════════════════
# MAIN
# ══════════════════════════════════════════════════════════════════════

def main():
    print("Phase 13: Reconciliation Tests")
    print("=" * 70)
    print("Do Phase 12 within-union findings contradict earlier papers?")
    print()

    tri = load_data()

    test1_results = test1_channel_switching(tri)
    test2_results = test2_si_regime_contingency(tri)
    write_summary(test1_results, test2_results)

    print("\n" + "=" * 70)
    print("Phase 13 complete.")
    print(f"Tables written to: {OUT_TABLES}")
    print("=" * 70)


if __name__ == '__main__':
    main()
